set more off

estimates clear
use heroin_deaths_state_month.dta, clear

quietly summ date if year==2010 & month==8
local reform = r(mean)

** quick aside to create a dataset of demographics
preserve
gen agege50_pop=age5064_pop+age6579_pop+agege80_pop
keep state fips date year month age2034_pop age3549_pop agege50_pop black* oth* hisp* urate 
gen quarter = 1 if month<=3
replace quarter = 2 if inrange(month,4,6)
replace quarter = 3 if inrange(month,7,9)
replace quarter = 4 if inrange(month,10,12)
	
collapse (mean) age2034_pop age3549_pop agege50_pop black* oth* hisp* urate ///
, by(state fips year quarter)
	
sort state year quarter
save ../factors/qtr_demographics.dta, replace
restore

	
** create variables that denote whether a state is in the higher or lower 
** risk group for a number of factors.
do risk_factors_1.do
	
gen trend = date - `reform'
gen treat = (date - `reform')*(date >= `reform')
replace trend=0 if date>=`reform'
   
gen index=trend

replace trend=0 if treat>0

xi i.month
gen hd_hs=rfact7_ge50==1 & rfact1_ge50==1
gen hd_ls=rfact7_ge50==1 & rfact1_ge50==0
gen ld_hs=rfact7_ge50==0 & rfact1_ge50==1
gen ld_ls=rfact7_ge50==0 & rfact1_ge50==0

sum hd_hs hd_ls ld_hs ld_ls

gen trend_hdhs=trend*hd_hs
gen trend_hdls=trend*hd_ls
gen trend_ldhs=trend*ld_hs
gen trend_ldls=trend*ld_ls

gen treat_hdhs=treat*hd_hs
gen treat_hdls=treat*hd_ls
gen treat_ldhs=treat*ld_hs
gen treat_ldls=treat*ld_ls

gen heroin_deaths_pc=heroin_deaths*100000/population
gen opioid_deaths_pc=opioid_deaths*100000/population
gen h_or_o_deaths_pc=h_or_o_deaths*100000/population

sum heroin_deaths_pc opioid_deaths_pc h_or_o_deaths_pc [aweight=population] if index<=-1 & index>=-12

gen agege50_pop=age5064_pop+age6579_pop+agege80_pop

	
** aside to make scatter plots of risk factors and outcomes in pre-period
preserve
keep if year>=2008 & year<2010
collapse (sum) heroin_deaths opioid_deaths h_or_o_deaths population ///
(max) orate_rf1 hrate_rf1, by(fips)

** convert from grams per 100k people to 1k people
replace orate_rf1 = orate_rf1/100

gen hrate = heroin_deaths/(population/100000)
gen orate = opioid_deaths/(population/100000)
gen horate = h_or_o_deaths/(population/100000)

** scale up heroin death rate risk factor to be per 100k people
replace hrate_rf1= 100000*hrate_rf1

reg orate orate_rf1, vce(robust)
reg orate orate_rf1 [aweight=population], vce(robust)

keep fips hrate-horate orate_rf1 hrate_rf1 
export delimited using ./scatterplot_maker.txt, replace
restore




	** Regressions
	** heroin
	areg heroin_deaths_pc trend_ldls treat_ldls trend_hdls treat_hdls trend_ldhs treat_ldhs trend_hdhs treat_hdhs ///
	_I* age2034_pop age3549_pop agege50_pop black* oth* hisp* urate ///
	[aweight=population] if year>=2004, abs(fips) vce(cluster fips)
	lincom treat_hdhs-trend_hdhs
	lincom treat_ldhs-trend_ldhs
	lincom treat_hdls-trend_hdls
	lincom treat_ldls-trend_ldls
	lincom treat_hdhs-trend_hdhs-(treat_ldls-trend_ldls)
	lincom treat_ldhs-trend_ldhs-(treat_ldls-trend_ldls)
	lincom treat_hdls-trend_hdls-(treat_ldls-trend_ldls)

	areg heroin_deaths_pc trend_ldls treat_ldls trend_hdls treat_hdls trend_ldhs treat_ldhs trend_hdhs treat_hdhs ///
	_I* age2034_pop age3549_pop agege50_pop black* oth* hisp* urate ///
	[aweight=population] if year>=2004 & year<=2012, abs(fips) vce(cluster fips)
	lincom treat_hdhs-trend_hdhs
	lincom treat_ldhs-trend_ldhs
	lincom treat_hdls-trend_hdls
	lincom treat_ldls-trend_ldls
	lincom treat_hdhs-trend_hdhs-(treat_ldls-trend_ldls)
	lincom treat_ldhs-trend_ldhs-(treat_ldls-trend_ldls)
	lincom treat_hdls-trend_hdls-(treat_ldls-trend_ldls)

	*** opioids
	areg opioid_deaths_pc trend_ldls treat_ldls trend_hdls treat_hdls trend_ldhs treat_ldhs trend_hdhs treat_hdhs ///
	_I* age2034_pop age3549_pop agege50_pop black* oth* hisp* urate ///
	[aweight=population] if year>=2004, abs(fips) vce(cluster fips)
	lincom treat_hdhs-trend_hdhs
	lincom treat_ldhs-trend_ldhs
	lincom treat_hdls-trend_hdls
	lincom treat_ldls-trend_ldls
	lincom treat_hdhs-trend_hdhs-(treat_ldls-trend_ldls)
	lincom treat_ldhs-trend_ldhs-(treat_ldls-trend_ldls)
	lincom treat_hdls-trend_hdls-(treat_ldls-trend_ldls)

	areg opioid_deaths_pc trend_ldls treat_ldls trend_hdls treat_hdls trend_ldhs treat_ldhs trend_hdhs treat_hdhs ///
	_I* age2034_pop age3549_pop agege50_pop black* oth* hisp* urate ///
	[aweight=population] if year>=2004 & year<=2012, abs(fips) vce(cluster fips)
	lincom treat_hdhs-trend_hdhs
	lincom treat_ldhs-trend_ldhs
	lincom treat_hdls-trend_hdls
	lincom treat_ldls-trend_ldls
	lincom treat_hdhs-trend_hdhs-(treat_ldls-trend_ldls)
	lincom treat_ldhs-trend_ldhs-(treat_ldls-trend_ldls)
	lincom treat_hdls-trend_hdls-(treat_ldls-trend_ldls)

	*** combined
	areg h_or_o_deaths_pc trend_ldls treat_ldls trend_hdls treat_hdls trend_ldhs treat_ldhs trend_hdhs treat_hdhs ///
	_I* age2034_pop age3549_pop agege50_pop black* oth* hisp* urate ///
	[aweight=population] if year>=2004, abs(fips) vce(cluster fips)
	lincom treat_hdhs-trend_hdhs
	lincom treat_ldhs-trend_ldhs
	lincom treat_hdls-trend_hdls
	lincom treat_ldls-trend_ldls
	lincom treat_hdhs-trend_hdhs-(treat_ldls-trend_ldls)
	lincom treat_ldhs-trend_ldhs-(treat_ldls-trend_ldls)
	lincom treat_hdls-trend_hdls-(treat_ldls-trend_ldls)

	areg h_or_o_deaths_pc trend_ldls treat_ldls trend_hdls treat_hdls trend_ldhs treat_ldhs trend_hdhs treat_hdhs ///
	_I* age2034_pop age3549_pop agege50_pop black* oth* hisp* urate ///
	[aweight=population] if year>=2004 & year<=2012, abs(fips) vce(cluster fips)
	lincom treat_hdhs-trend_hdhs
	lincom treat_ldhs-trend_ldhs
	lincom treat_hdls-trend_hdls
	lincom treat_ldls-trend_ldls
	lincom treat_hdhs-trend_hdhs-(treat_ldls-trend_ldls)
	lincom treat_ldhs-trend_ldhs-(treat_ldls-trend_ldls)
	lincom treat_hdls-trend_hdls-(treat_ldls-trend_ldls)

	summ heroin_deaths_pc opioid_deaths_pc h_or_o_deaths_pc ///
	[aweight=population] ///
	if (year==2009 & month>=8) | (year==2010 & month<8)
	
		
	** export data to excel for a couple graphs (hi-low graphs for each factor individuall)
	preserve
	collapse (sum) heroin_deaths opioid_deaths h_or_o_deaths population, ///
	by(year month date rfact1_ge50)
	gen h_drate = heroin_deaths/(population/100000)
	gen o_drate = opioid_deaths/(population/100000)
	gen h_or_o_drate = h_or_o_deaths/(population/100000)
	drop heroin_deaths opioid_deaths h_or_o_deaths population
	reshape wide h_drate o_drate h_or_o_drate, i(date year month) j(rfact1_ge50)
	foreach x in "h_drate" "o_drate" "h_or_o_drate" {
		rename `x'0 `x'_below
		rename `x'1 `x'_above
	}
	export delimited ./heroin_factor_graphing_data.txt, replace 
	restore
	
	preserve
	collapse (sum) heroin_deaths opioid_deaths h_or_o_deaths population, ///
	by(year month date rfact7_ge50)
	gen h_drate = heroin_deaths/(population/100000)
	gen o_drate = opioid_deaths/(population/100000)
	gen h_or_o_drate = h_or_o_deaths/(population/100000)
	drop heroin_deaths opioid_deaths h_or_o_deaths population
	reshape wide h_drate o_drate h_or_o_drate, i(date year month) j(rfact7_ge50)
	foreach x in "h_drate" "o_drate" "h_or_o_drate" {
		rename `x'0 `x'_below
		rename `x'1 `x'_above
	}
	export delimited ./opioid_factor_graphing_data.txt, replace 
	restore
	
	** and now export things for making the four group graph
	collapse (sum) heroin_deaths opioid_deaths h_or_o_deaths population, ///
	by(year month date rfact1_ge50 rfact7_ge50)
	drop if year<2004
	gen grp = 4 if rfact7_ge50==1 & rfact1_ge50==1
	replace grp= 2 if rfact7_ge50==1 & rfact1_ge50==0
	replace grp= 3 if rfact7_ge50==0 & rfact1_ge50==1
	replace grp= 1 if rfact7_ge50==0 & rfact1_ge50==0
	drop rfact1_ge50 rfact7_ge50
	gen hrate = heroin_deaths/(population/100000)
	drop if heroin_deaths<10
	drop heroin_deaths opioid_deaths h_or_o population 
	reshape wide hrate, i(date year month) j(grp)
	rename hrate1 hrate_ll
	rename hrate2 hrate_hl
	rename hrate3 hrate_lh
	rename hrate4 hrate_hh
	export delimited ./four_group_graphiing_data.txt, replace 
	

	




		
